% This code can be executed to generate Figure 6. The purple curve is 
% theoretical peeling maximum (PeelF) is expressed as a function of the 
% junction angle. The yellow curve is the elastic force (intersectUFelas) 
% where the peeling and elastic energies are equal (intersectU). The orange 
% curve is point where the elastic model becomes linear at x=alpha 
% (linearpts). The tensile strength of the tape (TS_tape) is the black 
% dashed line. The lap shear strength (lapshearstrength) is plotted as a 
% solid black line. The pink curve is the average load for junctions with a
% G2 (R=2.6 inch) pre-buckle (datloc2). The blue-green curve is the average load for
% junctions with a G1 (R=4 inch) pre-buckle (datloc4). Additional formatting
% is done externally to match the published Firgure 6.
% The folder titled, "Raw Data Test Files," must be open to run this code!
clear;
%Tape Constants
t=0.058; %mm
w=12.7; %mm
E=467.9745; % Pa
gamma=10.3014/1000;  %N/mm
beta=mean([189.393 190.444 184.519 182.807 186.996]);
TS_tape=23.2416;
slope=(E*w)./1000; %N/mm

%Determine Our Model (Energy Equivalence and Transition to Linearity
%(alpha))
%initialize variables
angles=linspace(1.25,76,10000);
intersectUx=zeros(1,length(angles));
intersectUFelas=zeros(1,length(angles));
intersectU=zeros(1,length(angles));
intersectFx=zeros(1,length(angles));
xmaxs=zeros(1,length(angles));
TS_tapeFelasx=zeros(1,length(angles));
linearpts=zeros(1,length(angles));

for i=1:length(angles)
    phi=angles(i);
    xmax=2.*w.*tand(phi./2);
    xmaxs(i)=xmax;
    %Peeling Energy
    Upeelinc=@(x) 2.*gamma.*((cscd(phi./2)).^2).*(1./2).*(x.^2).*(cotd(phi./2));
    Upeeldec=@(x) 2.*gamma.*((cscd(phi./2)).^2).*((1)./2).*((-1.*(xmax-x).^2)+2.*(xmax./2).^2).*(cotd(phi./2));
    
    xx=linspace(0,20,10000);
    for j=1:length(xx)
        pos=xx(j);
        if pos>=(xmax/2)
            xxinc=xx(1:j);
            Upeelincvec=Upeelinc(xxinc);
            half=j;
            break
        end
    end
    for h=half:length(xx)
        pos=xx(h);
        if pos>=xmax
            xxdec=xx(half:h);
            Upeeldecvec=Upeeldec(xxdec);
            break
        end
    end
    %add tail of constant value to make sure vectors are same length
    tailU=ones(1,(length(xx)-h)-1).*max(Upeeldecvec);
    Upeel=[Upeelincvec Upeeldecvec tailU];
    
    %Elastic Force
    Felas=@(x) -(1/2).*E.*t.*x.*tand(phi).*log(abs((2.*x.*cotd(phi./2)-beta.*tand(phi))./(-beta.*tand(phi))));
    Felas_val=Felas(xx);
    for n=1:(length(Felas_val)-1)
        m=(Felas_val(n+1)-Felas_val(n))/(xx(n+1)-xx(n));
        b=Felas_val(n)-m*xx(n);
        %linearize elastic force
        if m>slope
            xremain=linspace(xx(n),max(xx),(length(xx)-n+1));
            Felas_val(n:length(xx))=xremain.*slope+b;
            break
        end
        linearpts(i)=xx(n);
    end
    %Elastic Energy
    Uelas=cumtrapz(xx,Felas_val);
    
    %Difference between elastic and peeling energies to determine where
    %energies intersect
    subU=Uelas-Upeel;
    
    for k=1:length(subU)-1
        test=subU(k);
        if test>0
            intersectUx(i)=xx(k);
            intersectUFelas(i)=Felas_val(k);
            intersectU(i)=Uelas(k);
            break
        end
    end
    
    Fpeelinc=@(x) 2.*gamma.*(x).*((cscd(phi./2)).^2).*(cotd(phi./2));
    Fpeeldec=@(x) 2.*gamma.*(xmax-(x)).*((cscd(phi./2)).^2).*(cotd(phi./2));
    
    tailF=zeros(1,length(tailU));
    Fpeel=[Fpeelinc(xxinc) Fpeeldec(xxdec) tailF];
    subF=Felas_val-Fpeel;
    
    %Determine Alpha
    linearptsFelas(i)=Felas_val(n);
    for p=1:length(subF)-1
        test=subF(p);
        if test>0
            intersectFx(i)=xx(p);
            break
        end
    end
    for d=1:length(Felas_val)
        test=Felas_val(d);
        if test>TS_tape
            TS_tapeFelasx(i)=xx(d);
            break
        end
    end
end
anglestheo=angles;

%Interperet Raw Data: different sets of data require different treatment
%due to the organization in the .txt files and category of the sample that
%the file cooresponds to.

%Tensile Tester Trials for 4" prebuckle
filenameangles4=[flip(linspace(20,90,15)),flip(linspace(22.5,87.5,14)),flip(linspace(20,90,15)),flip(linspace(22.5,87.5,14)), 27.5, 27.5];
trials4=linspace(109,168,60);
names4=[];
for l=1:29
    testname=['2021-10-26_',num2str(trials4(l)),'_',num2str(round(filenameangles4(l),1)),'.txt'];
    if str2num(testname(17))==0
        testname=['2021-10-26_',num2str(trials4(l)),'_',num2str(round(filenameangles4(l),1)),'.0.txt'];
    elseif str2num(testname(17))==5
        testname=['2021-10-26_',num2str(trials4(l)),'_',num2str(round(filenameangles4(l),1)),'.0.txt'];
    end
    names4=[names4 testname];
end
for l=30:length(trials4)
    testname=['2021-10-28_',num2str(trials4(l)),'_',num2str(round(filenameangles4(l),1)),'.txt'];
    if str2num(testname(17))==0
        testname=['2021-10-28_',num2str(trials4(l)),'_',num2str(round(filenameangles4(l),1)),'.0.txt'];
    elseif str2num(testname(17))==5
        testname=['2021-10-28_',num2str(trials4(l)),'_',num2str(round(filenameangles4(l),1)),'.0.txt'];
    end
    names4=[names4 testname];
end

%short algorithm to determine where each file name starts and ends to
%determine the junction angle information which is stored in the filenames
for i=1:length(names4)
    chunk=i+3;
    test=names4(i:chunk) == '.txt';
    if all(test)
        numchar=chunk;  % number of characters in each file name
        break
    end
end

%initialize variables
ntrials4=length(names4)/numchar;
jj4=zeros(1,ntrials4);
jj4(1)=1;
kk4=zeros(1,ntrials4);
kk4(1)=numchar;
elasticmax4=[];
elasticmaxangles4=[];
peelmax4=[];
peelmaxangles4=[];
peakload4=[];
angles4=[];
for i=1:ntrials4
    clear Data displacement load time
    filename=names4(jj4(i):kk4(i));
    fileID=importdata(filename);
    j=jj4(i)+numchar;
    k=kk4(i)+numchar;
    jj4(i+1)=j;
    kk4(i+1)=k;
    Data=fileID.data;
    displacement=Data(:,1);  %mm
    load=Data(:,2);  %N
    adjusttrials=[109 110 112 113 114 115 117 118 123 124 125 128 129 130 131 132 134 136 137 138 139 140 141 142 143 148 153 154 155 156 157 158 159 160 161 163];
    adjusttripletroubletrials=[81 112 113 134 140 141 143 153 157 160];
    adjustquadtroubletrials=[19];
    adjustbigtroubletrials=[];
    phi=str2num(filename(16:19));
    %Determine peaks of tensile data
    if any(ismember(adjusttrials,str2num(filename(12:14)))) %these trials have 2 peaks, global: elastic, local: peeling
        [pks, pkloc]=findpeaks(load,displacement);
        [pks2, pkloc2]=findpeaks(pks,pkloc);
        allpks=[pks2 pkloc2];
        sortpks=sortrows(allpks,'descend');
        peelmax4=[peelmax4 sortpks(2,1)];
        peelmaxangles4=[peelmaxangles4 phi];
        elasticmax4=[elasticmax4 sortpks(1,1)];
        elasticmaxangles4=[elasticmaxangles4 phi];
        if any(ismember(adjusttripletroubletrials,str2num(filename(12:14)))) 
            peelmax4(length(peelmax4))=sortpks(3,1);  %these trials have the issues discussed above
        end
        if any(ismember(adjustquadtroubletrials,str2num(filename(12:14))))
            peelmax4(length(peelmax4))=sortpks(4,1);
        end
        if any(ismember(adjustbigtroubletrials,str2num(filename(12:14))))
            [pks3, pkloc3]=findpeaks(pks2,pkloc2);
            allpks=[pks3 pkloc3];
            sortpks=sortrows(allpks,'descend');
            if str2num(filename(12:14))==2
                peelmax4(length(peelmax4))=sortpks(2,1);
            end
            if str2num(filename(12:14))==3
                peelmax4(length(peelmax4))=sortpks(3,1);
            end
        end
    else
        [pks, pkloc]=findpeaks(load,displacement);  %these trials are peeling only
        allpks=[pks pkloc];
        sortpks=sortrows(allpks,'descend');
        peakload4=[peakload4 sortpks(1,1)];
        angles4=[angles4 phi];
    end
end

% Analysis
%organize 4" prebuckle data into one matrix
dat4=[angles4; peakload4];
dat4=sortrows(dat4',1);
elasticmaxdat4=[elasticmaxangles4; elasticmax4];
peelmaxdat4=[peelmaxangles4; peelmax4];

%Determine Mean and Range for given sections of the peeled data to draw the 
%Mean and Range lines and regions respectively
rangedat4=sortrows([dat4(:,1)' elasticmaxdat4(1,:); dat4(:,2)'./areacalculation(w,dat4(:,1)') elasticmaxdat4(2,:)./areacalculation(w,elasticmaxdat4(1,:))]',1);

%binwidth=(max(rangedat(:,1))-min(rangedat(:,1)))/numbin;
%binlims=linspace(min(rangedat(:,1)),max(rangedat(:,1)),numbin);
binlims4=[linspace(min(rangedat4(:,1)),20,5) linspace(22,30,4) linspace(33,max(rangedat4(:,1)),9)];
numbin4=length(binlims4);
datloc4=zeros(6,numbin4-1);
for i=1:(numbin4-1)
    minlim=binlims4(i);
    maxlim=binlims4(i+1);
    if i==1
        datloc4(1,i)=minlim;
    elseif i==(numbin4-1)
        datloc4(1,i)=maxlim;
    else
        datloc4(1,i)=(maxlim+minlim)/2;
    end
    bin=[];
    for j=1:length(rangedat4(:,1))
        angle=rangedat4(j,1);
        if angle>=minlim && angle<=maxlim
            bin=[bin rangedat4(j,2)];
        end
    end
    datloc4(2,i)=mean(bin);
    datloc4(3,i)=std(bin);
    if length(bin)==1
        datloc4(4,i)=mean(bin);
        datloc4(5,i)=mean(bin);
    else
        datloc4(4,i)=min(bin);
        datloc4(5,i)=max(bin);
    end
    datloc4(6,i)=length(bin);
    %datloc(7,i)=datloc(2,i)-datloc(3,i);
    %datloc(8,i)=datloc(2,i)+datloc(3,i);
end

%Plotting
figure;
hold on

%Plot Average for 4" Data Points
plot(datloc4(1,:),datloc4(2,:),'Marker','none','Color',[0/220 150/255 150/255],'LineStyle','-','LineWidth',3);

%Optional: plot range for 4" Data points
%jbfill(datloc4(1,:),datloc4(4,:),datloc4(5,:),[106/220 255/255 255/255],[1 1 1],[0 0 0],119/255);
%plot(datloc4(1,:),datloc4(4,:),'Marker','none','Color',[0/255 150/255 150/255],'LineStyle','-');
%plot(datloc4(1,:),datloc4(5,:),'Marker','none','Color',[0/255 150/255 150/255],'LineStyle','-');

%split 4" prebuckle data into nestled (n) and locked (l) data sets
nlocked4=0;
nnestled4=0;
R3loads4=[];
R2loads4=[];
%Plot Data Points based on conformation and failure regions
%optional: scatter plot of all data points
for c=1:length(dat4(:,1))
    phi=dat4(c,1);
    if phi<21.9
        %scatter(dat4(c,1),(dat4(c,2)./areacalculation(w,dat4(c,1))'),30,'filled','MarkerFaceAlpha',1/3,'MarkerFaceColor',[0 0 0]...
        %    ,'MarkerEdgeColor',[0 0 0],'MarkerEdgeAlpha',2/3);
    elseif phi<=30 && dat4(c,2)/areacalculation(w,dat4(c,1))>=0.0256
        %scatter(dat4(c,1),(dat4(c,2)./areacalculation(w,dat4(c,1))'),30,'filled','MarkerFaceAlpha',1/3,'MarkerFaceColor',[255/255 112/255 101/255]...
        %    ,'MarkerEdgeColor',[0 0 0],'MarkerEdgeAlpha',2/3);
        nlocked4=nlocked4+1;
        R2loads4=[R2loads4 dat4(c,2)./areacalculation(w,dat4(c,1))];
    elseif phi<=30 && dat4(c,2)/areacalculation(w,dat4(c,1))<0.0256
        %scatter(dat4(c,1),(dat4(c,2)./areacalculation(w,dat4(c,1))'),30,'filled','MarkerFaceAlpha',1/3,'MarkerFaceColor',[0/255 190/255 30/255]...
        %    ,'MarkerEdgeColor',[0 0 0],'MarkerEdgeAlpha',2/3);
        nnestled4=nnestled4+1;
        R2loads4=[R2loads4 dat4(c,2)./areacalculation(w,dat4(c,1))];
    elseif phi>30 && dat4(c,2)/areacalculation(w,dat4(c,1))<0.037
        %scatter(dat4(c,1),(dat4(c,2)./areacalculation(w,dat4(c,1))'),30,'filled','MarkerFaceAlpha',1/3,'MarkerFaceColor',[0/255 190/255 30/255]...
        %    ,'MarkerEdgeColor',[0/255 100/255 0/255],'MarkerEdgeAlpha',2/3);
        nnestled4=nnestled4+1;
        R3loads4=[R3loads4 dat4(c,2)./areacalculation(w,dat4(c,1))];
    else
        %scatter(dat4(c,1),(dat4(c,2)./areacalculation(w,dat4(c,1))'),30,'filled','MarkerFaceAlpha',128/255,'MarkerFaceColor',[255/255 112/255 101/255]...
        %    ,'MarkerEdgeColor',[180/255 0/255 0/255],'MarkerEdgeAlpha',128/255);
        nlocked4=nlocked4+1;
        R3loads4=[R3loads4 dat4(c,2)./areacalculation(w,dat4(c,1))];
    end
end
%scatter(elasticmaxdat4(1,:)',elasticmaxdat4(2,:)'./areacalculation(w,elasticmaxdat4(1,:))',30,'filled','MarkerFaceAlpha',128/255,'MarkerFaceColor',[255/255 112/255 101/255]...
%    ,'MarkerEdgeColor',[180/255 0/255 0/255],'MarkerEdgeAlpha',128/255);

R3loads4=[R3loads4 elasticmaxdat4(2,:)./areacalculation(w,elasticmaxdat4(1,:))];

%scatter(peelmaxdat4(1,:)',peelmaxdat4(2,:)'./areacalculation(w,peelmaxdat4(1,:))',30,'filled','MarkerFaceAlpha',0/3,'MarkerFaceColor',[180/255 0/255 0/255]...
%    ,'MarkerEdgeColor',[180/255 0/255 0/255],'MarkerEdgeAlpha',3/3);

%determine average strengths for nestled and locked data sets for the 4"
%prebuckled data
nlocked4=nlocked4+length(elasticmaxdat4);
percentlocked4=nlocked4/(nlocked4+nnestled4);
R3avgload4=mean(R3loads4);
R2avgload4=mean(R2loads4);
R23avgload4=mean([R3loads4 R2loads4]);


%Tensile Tester Trials for 2.5" preload
filenameangles2=[90,flip(linspace(67.5,85,8)),flip(linspace(20,62.5,18)),65,87.5,85,82.5,77.5,flip(linspace(20,70,21)),87.5,75,72.5];
%leaving out test runs: 169, 171, 180, 201, 202, 205, 207, 208, 230, 232
%these trials were not executed correctly in the tensile tester
trials2=[170,linspace(172,179,8),linspace(181,200,20),203,204,206,linspace(209,229,21),231,233,234];
names2=[];
for l=1:length(trials2)
    trial=trials2(l);
    testname=['2021-11-03_',num2str(trial),'_',num2str(round(filenameangles2(l),1)),'.txt'];
    if str2num(testname(17))==0
        testname=['2021-11-03_',num2str(trial),'_',num2str(round(filenameangles2(l),1)),'.0.txt'];
    elseif str2num(testname(17))==5
        testname=['2021-11-03_',num2str(trial),'_',num2str(round(filenameangles2(l),1)),'.0.txt'];
    end
    names2=[names2 testname];
end

%short algorithm to determine where each file name starts and ends to
%determine the junction angle information which is stored in the filenames
for i=1:length(names2)
    chunk=i+3;
    test=names2(i:chunk) == '.txt';
    if all(test)
        numchar=chunk;  % number of characters in each file name
        break
    end
end
%initialize variables
ntrials2=length(names2)/numchar;
jj2=zeros(1,ntrials2);
jj2(1)=1;
kk2=zeros(1,ntrials2);
kk2(1)=numchar;
elasticmax2=[];
elasticmaxangles2=[];
peelmax2=[];
peelmaxangles2=[];
peakload2=[];
angles2=[];
for i=1:ntrials2
    clear Data displacement load time
    filename=names2(jj2(i):kk2(i));
    fileID=importdata(filename);
    j=jj2(i)+numchar;
    k=kk2(i)+numchar;
    jj2(i+1)=j;
    kk2(i+1)=k;
    Data=fileID.data;
    displacement=Data(:,1);  %mm
    load=Data(:,2);  %N
    adjusttrials=[170 172 173 174 175 176 177 178 179 181 182 183 184 185 186 187 188 189 190 191 199 200 203 204 206 209 210 211 212 213 214 215 216 217 218 219 220 221 231 233 234];
    adjusttripletroubletrials=[176 182 188 189 190 199 204 209 210 211 215 217 219];
    adjustquadtroubletrials=[19];
    adjustbigtroubletrials=[];
    phi=str2num(filename(16:19));
    %determine peaks for 2.5" prebuckle data
    if any(ismember(adjusttrials,str2num(filename(12:14)))) %these trials have 2 peaks, global: elastic, local: peeling
        [pks, pkloc]=findpeaks(load,displacement);
        [pks2, pkloc2]=findpeaks(pks,pkloc);
        allpks=[pks2 pkloc2];
        sortpks=sortrows(allpks,'descend');
        peelmax2=[peelmax2 sortpks(2,1)];
        peelmaxangles2=[peelmaxangles2 phi];
        elasticmax2=[elasticmax2 sortpks(1,1)];
        elasticmaxangles2=[elasticmaxangles2 phi];
        if any(ismember(adjusttripletroubletrials,str2num(filename(12:14)))) 
            peelmax2(length(peelmax2))=sortpks(3,1);  %these trials have the issues discussed above
        end
        if any(ismember(adjustquadtroubletrials,str2num(filename(12:14))))
            peelmax2(length(peelmax2))=sortpks(4,1);
        end
        if any(ismember(adjustbigtroubletrials,str2num(filename(12:14))))
            [pks3, pkloc3]=findpeaks(pks2,pkloc2);
            allpks=[pks3 pkloc3];
            sortpks=sortrows(allpks,'descend');
            if str2num(filename(12:14))==2
                peelmax2(length(peelmax2))=sortpks(2,1);
            end
            if str2num(filename(12:14))==3
                peelmax2(length(peelmax2))=sortpks(3,1);
            end
        end
    else
        [pks, pkloc]=findpeaks(load,displacement);  %these trials are peeling only
        allpks=[pks pkloc];
        sortpks=sortrows(allpks,'descend');
        peakload2=[peakload2 sortpks(1,1)];
        angles2=[angles2 phi];
    end
end

% Analysis
%organize 2.5" prebuckle data into one matrix
dat2=[angles2; peakload2];
dat2=sortrows(dat2',1);
elasticmaxdat2=[elasticmaxangles2; elasticmax2];
peelmaxdat2=[peelmaxangles2; peelmax2];

%Determine Mean and Range for given sections of the peeled data to draw the 
%Mean and Range lines and regions respectively
rangedat2=sortrows([dat2(:,1)' elasticmaxdat2(1,:); dat2(:,2)'./areacalculation(w,dat2(:,1)') elasticmaxdat2(2,:)./areacalculation(w,elasticmaxdat2(1,:))]',1);

%binwidth=(max(rangedat(:,1))-min(rangedat(:,1)))/numbin;
%binlims=linspace(min(rangedat(:,1)),max(rangedat(:,1)),numbin);
binlims2=[linspace(min(rangedat2(:,1)),20,5) linspace(22,30,4) linspace(33,max(rangedat2(:,1)),9)];
numbin2=length(binlims2);
datloc2=zeros(6,numbin2-1);
for i=1:(numbin2-1)
    minlim=binlims2(i);
    maxlim=binlims2(i+1);
    if i==1
        datloc2(1,i)=minlim;
    elseif i==(numbin2-1)
        datloc2(1,i)=maxlim;
    else
        datloc2(1,i)=(maxlim+minlim)/2;
    end
    bin=[];
    for j=1:length(rangedat2(:,1))
        angle2=rangedat2(j,1);
        if angle2>=minlim && angle2<=maxlim
            bin=[bin rangedat2(j,2)];
        end
    end
    datloc2(2,i)=mean(bin);
    datloc2(3,i)=std(bin);
    if length(bin)==1
        datloc2(4,i)=mean(bin);
        datloc2(5,i)=mean(bin);
    else
        datloc2(4,i)=min(bin);
        datloc2(5,i)=max(bin);
    end
    datloc2(6,i)=length(bin);
    %datloc(7,i)=datloc(2,i)-datloc(3,i);
    %datloc(8,i)=datloc(2,i)+datloc(3,i);
end

%Plot Range and Average for 2.5" Data Points
plot(datloc2(1,:),datloc2(2,:),'Marker','none','Color',[255/255 150/255 150/255],'LineStyle','-','LineWidth',3);

%Optional: plot range for 2.5" data points
%jbfill(datloc2(1,:),datloc2(4,:),datloc2(5,:),[106/220 255/255 255/255],[1 1 1],[0 0 0],119/255);
%plot(datloc2(1,:),datloc2(4,:),'Marker','none','Color',[0/255 150/255 150/255],'LineStyle','-');
%plot(datloc2(1,:),datloc2(5,:),'Marker','none','Color',[0/255 150/255 150/255],'LineStyle','-');

nlocked2=0;
nnestled2=0;
R3loads2=[];
R2loads2=[];
%Plot Data Points based on conformation and failure regions
%optional: scatter plot all data points
for c=1:length(dat2(:,1))
    phi=dat2(c,1);
    if phi<21.9
        %scatter(dat2(c,1),(dat2(c,2)./areacalculation(w,dat2(c,1))'),30,'filled','MarkerFaceAlpha',1/3,'MarkerFaceColor',[0 0 0]...
        %    ,'MarkerEdgeColor',[0 0 0],'MarkerEdgeAlpha',2/3);
    elseif phi<=30 && dat2(c,2)/areacalculation(w,dat2(c,1))>=0.0256
        %scatter(dat2(c,1),(dat2(c,2)./areacalculation(w,dat2(c,1))'),30,'filled','MarkerFaceAlpha',1/3,'MarkerFaceColor',[255/255 112/255 101/255]...
        %    ,'MarkerEdgeColor',[0 0 0],'MarkerEdgeAlpha',2/3);
        nlocked2=nlocked2+1;
        R2loads2=[R2loads2 dat2(c,2)./areacalculation(w,dat2(c,1))];
    elseif phi<=30 && dat2(c,2)/areacalculation(w,dat2(c,1))<0.0256
        %scatter(dat2(c,1),(dat2(c,2)./areacalculation(w,dat2(c,1))'),30,'filled','MarkerFaceAlpha',1/3,'MarkerFaceColor',[0/255 190/255 30/255]...
        %    ,'MarkerEdgeColor',[0 0 0],'MarkerEdgeAlpha',2/3);
        nnestled2=nnestled2+1;
        R2loads2=[R2loads2 dat2(c,2)./areacalculation(w,dat2(c,1))];
    elseif phi>30 && dat2(c,2)/areacalculation(w,dat2(c,1))<0.037
        %scatter(dat2(c,1),(dat2(c,2)./areacalculation(w,dat2(c,1))'),30,'filled','MarkerFaceAlpha',1/3,'MarkerFaceColor',[0/255 190/255 30/255]...
        %    ,'MarkerEdgeColor',[0/255 100/255 0/255],'MarkerEdgeAlpha',2/3);
        nnestled2=nnestled2+1;
        R3loads2=[R3loads2 dat2(c,2)./areacalculation(w,dat2(c,1))];
    else
        %scatter(dat2(c,1),(dat2(c,2)./areacalculation(w,dat2(c,1))'),30,'filled','MarkerFaceAlpha',128/255,'MarkerFaceColor',[255/255 112/255 101/255]...
        %    ,'MarkerEdgeColor',[180/255 0/255 0/255],'MarkerEdgeAlpha',128/255);
        nlocked2=nlocked2+1;
        R3loads2=[R3loads2 dat2(c,2)./areacalculation(w,dat2(c,1))];
    end
end
%scatter(elasticmaxdat2(1,:)',elasticmaxdat2(2,:)'./areacalculation(w,elasticmaxdat2(1,:))',30,'filled','MarkerFaceAlpha',128/255,'MarkerFaceColor',[255/255 112/255 101/255]...
%    ,'MarkerEdgeColor',[180/255 0/255 0/255],'MarkerEdgeAlpha',128/255);

R3loads2=[R3loads2 elasticmaxdat2(2,:)./areacalculation(w,elasticmaxdat2(1,:))];

%scatter(peelmaxdat2(1,:)',peelmaxdat2(2,:)'./areacalculation(w,peelmaxdat2(1,:))',30,'filled','MarkerFaceAlpha',0/3,'MarkerFaceColor',[180/255 0/255 0/255]...
%    ,'MarkerEdgeColor',[180/255 0/255 0/255],'MarkerEdgeAlpha',3/3);

%determine average strengths for 2.5" prebuckled data 
nlocked2=nlocked2+length(elasticmaxdat2);
percentlocked2=nlocked2/(nlocked2+nnestled2);
R3avgload2=mean(R3loads2);
R2avgload2=mean(R2loads2);
R23avgload2=mean([R3loads2 R2loads2]);

%Tensile Tester Trials of non-prebuckled data (same data as Figure 5A)
names=['2019-02-08_001_70.0.txt',...
    '2019-02-27_002_75.0.txt',...
    '2019-02-27_003_66.0.txt',...
    '2019-02-27_004_65.0.txt',...
    '2019-02-27_005_52.6.txt',...
    '2019-02-27_006_40.9.txt',...
    '2019-02-27_007_28.0.txt',...
    '2019-03-01_008_85.1.txt',...
    '2019-03-01_009_52.0.txt',...
    '2019-03-01_010_49.5.txt',...
    '2019-03-01_011_30.6.txt',...
    '2020-04-15_001_10.0.txt',...
    '2020-04-15_002_12.5.txt',...
    '2020-04-15_003_15.0.txt',...
    '2020-04-15_004_17.5.txt',...
    '2020-04-15_005_20.0.txt',...
    '2020-04-15_006_22.5.txt',...
    '2020-04-15_007_25.0.txt',...
    '2020-04-15_008_27.5.txt',...
    '2020-04-15_009_30.0.txt',...
    '2020-04-15_010_32.5.txt',...
    '2020-04-15_011_35.0.txt',...
    '2020-04-15_012_37.5.txt',...
    '2020-04-15_013_40.0.txt',...
    '2020-04-15_014_42.5.txt',...
    '2020-04-15_015_45.0.txt',...
    '2020-04-15_016_47.5.txt',...
    '2020-04-15_017_50.0.txt',...
    '2020-04-15_018_52.5.txt',...
    '2020-04-15_019_55.0.txt',...
    '2020-04-15_020_57.5.txt',...
    '2020-04-15_021_60.0.txt',...
    '2020-04-15_022_65.0.txt',...
    '2020-04-15_023_70.0.txt',...
    '2020-04-15_024_75.0.txt',...
    '2020-05-04_025_10.0.txt',...
    '2020-05-04_026_12.5.txt',...
    '2020-05-04_027_15.0.txt',...
    '2020-05-04_028_17.5.txt',...
    '2020-05-04_029_20.0.txt',...
    '2020-05-04_030_22.5.txt',...
    '2020-05-04_031_25.0.txt',...
    '2020-05-04_032_27.5.txt',...
    '2020-05-04_033_30.0.txt',...
    '2020-05-04_034_17.5.txt',...
    '2020-05-04_035_20.0.txt',...
    '2020-05-04_036_22.5.txt',...
    '2020-05-04_037_25.0.txt',...
    '2020-05-04_038_27.5.txt',...
    '2020-05-04_039_30.0.txt',...
    '2020-05-04_040_20.0.txt',...
    '2020-05-04_041_22.5.txt',...
    '2020-05-04_042_25.0.txt',...
    '2020-05-04_043_27.5.txt',... 
    '2020-05-12_046_17.5.txt',...
    '2020-05-12_047_20.0.txt',...
    '2020-05-12_048_22.5.txt',...
    '2020-05-12_049_25.0.txt',...
    '2020-05-12_050_27.5.txt',...
    '2020-05-12_051_30.0.txt',...
    '2020-05-12_052_17.5.txt',...
    '2020-05-12_053_20.0.txt',...
    '2020-05-12_054_22.5.txt',...
    '2020-05-12_055_25.0.txt',...
    '2020-05-12_056_27.5.txt',...
    '2020-05-12_057_30.0.txt',...
    '2020-05-12_058_17.5.txt',...
    '2020-05-12_059_20.0.txt',...
    '2020-05-12_060_22.5.txt',...
    '2020-05-12_061_25.0.txt',...
    '2020-05-12_062_27.5.txt',...
    '2020-05-12_063_30.0.txt',...
    '2020-05-12_064_17.5.txt',...
    '2020-05-12_065_20.0.txt',...
    '2020-05-12_066_22.5.txt',...
    '2020-05-12_067_25.0.txt',...
    '2020-05-12_068_27.5.txt',...
    '2020-05-12_069_30.0.txt',...
    '2020-05-14_070_10.0.txt',...
    '2020-05-14_071_12.5.txt',...
    '2020-05-14_072_15.0.txt',...
    '2020-05-14_073_10.0.txt',...
    '2020-05-14_074_12.5.txt',...
    '2020-05-14_075_15.0.txt',...
    '2020-05-14_076_17.5.txt',...
    '2020-05-21_077_15.0.txt',...
    '2020-05-21_078_42.5.txt',...
    '2020-05-21_079_35.0.txt',...
    '2020-05-21_080_50.0.txt',...
    '2020-05-21_081_37.5.txt',...
    '2020-05-21_082_45.0.txt',...
    '2020-05-21_083_47.5.txt',...
    '2020-05-21_084_32.5.txt',...
    '2020-05-21_085_40.0.txt',...
    '2020-05-21_086_22.5.txt',...
    '2021-03-16_100_50.0.txt',...
    '2021-03-16_101_15.0.txt',...
    '2021-03-16_102_15.0.txt',...
    '2021-03-16_103_15.0.txt',...
    '2021-03-16_104_17.5.txt',...
    '2021-03-16_105_10.0.txt'];

%short algorithm to determine where each file name starts and ends to
%determine the junction angle information which is stored in the filenames
for i=1:length(names)
    chunk=i+3;
    test=names(i:chunk) == '.txt';
    if all(test)
        numchar=chunk;  % number of characters in each file name
        break
    end
end
%initialize variables
ntrials=length(names)/numchar;
jj=zeros(1,ntrials);
jj(1)=1;
kk=zeros(1,ntrials);
kk(1)=numchar;
elasticmax=[];
elasticmaxangles=[];
peelmax=[];
peelmaxangles=[];
peakload=[];
angles=[];
for i=1:ntrials
    clear Data displacement load time
    filename=names(jj(i):kk(i));
    fileID=importdata(filename);
    j=jj(i)+numchar;
    k=kk(i)+numchar;
    jj(i+1)=j;
    kk(i+1)=k;
    Data=fileID.data;
    if str2num(filename(4))==9  %older data sets require different treatment
        displacement=Data(:,1).*(50E-3);  %mm
        load=abs(Data(:,2)-abs(max(Data(:,2))));  %N
        adjusttrials=[2 3 5 6 10];       %we need to pick out the elastic and peeling peaks.
        adjusttripletroubletrials=[10];  %theoretically the highest peak is the elastic peak and the second is the peeling
        adjustquadtroubletrials=[];      %however, some trials have multiple peaks in the data so we basically manually pickout
        adjustbigtroubletrials=[2 3];    %the proper peaks. Only 5 samples have this issue and are labeled by the 
    else                                 %"adjust#troubletrials" varialbes
        displacement=Data(:,1);  %mm
        load=Data(:,2);  %N
        adjusttrials=[13 15 17 18 19 20 22 24 81 82 85];
        adjusttripletroubletrials=[81];
        adjustquadtroubletrials=[19];
        adjustbigtroubletrials=[];
    end
    phi=str2num(filename(16:19));
    %determine peaks for non-prebuckled data
    if any(ismember(adjusttrials,str2num(filename(12:14)))) %these trials have 2 peaks, global: elastic, local: peeling
        [pks, pkloc]=findpeaks(load,displacement);
        [pks2, pkloc2]=findpeaks(pks,pkloc);
        allpks=[pks2 pkloc2];
        sortpks=sortrows(allpks,'descend');
        peelmax=[peelmax sortpks(2,1)];
        peelmaxangles=[peelmaxangles phi];
        elasticmax=[elasticmax sortpks(1,1)];
        elasticmaxangles=[elasticmaxangles phi];
        if any(ismember(adjusttripletroubletrials,str2num(filename(12:14)))) 
            peelmax(length(peelmax))=sortpks(3,1);  %these trials have the issues discussed above
        end
        if any(ismember(adjustquadtroubletrials,str2num(filename(12:14))))
            peelmax(length(peelmax))=sortpks(4,1);
        end
        if any(ismember(adjustbigtroubletrials,str2num(filename(12:14))))
            [pks3, pkloc3]=findpeaks(pks2,pkloc2);
            allpks=[pks3 pkloc3];
            sortpks=sortrows(allpks,'descend');
            if str2num(filename(12:14))==2
                peelmax(length(peelmax))=sortpks(2,1);
            end
            if str2num(filename(12:14))==3
                peelmax(length(peelmax))=sortpks(3,1);
            end
        end
    else
        [pks, pkloc]=findpeaks(load,displacement);  %these trials are peeling only
        allpks=[pks pkloc];
        sortpks=sortrows(allpks,'descend');
        peakload=[peakload sortpks(1,1)];
        angles=[angles phi];
    end
end

%Peak Load Only Trials
peakloaddata=importdata('Peak Load Only Raw Data.txt');
peakloaddata2=peakloaddata.data;
anglesexpinv=peakloaddata2(:,1);
maxloadexp=flip(peakloaddata2(:,2))';
anglesexp=flip(anglesexpinv)';

% Analysis
%organize non-prebuckled data into one matrix
dat=[anglesexp angles; maxloadexp./areacalculation(w,anglesexp) peakload./areacalculation(w,angles)];
dat=sortrows(dat',1);
elasticmaxdat=[elasticmaxangles; elasticmax./areacalculation(w,elasticmaxangles)];
peelmaxdat=[peelmaxangles; peelmax./areacalculation(w,peelmaxangles)];

%Peeling Model
PeelF=@(phi) 2.*gamma.*w.*((cscd(phi./2)).^2)+0.4536;
phiphi=linspace(min(dat(:,1)),90,10000);
peel_f=PeelF(anglestheo)./areacalculation(w,anglestheo);

%Determine Boundary Between Regions II and III and Regions I and II
linearptsFelas=linearptsFelas./areacalculation(w,anglestheo);
intersectUFelas=intersectUFelas./areacalculation(w,anglestheo);
for i=1:length(anglestheo)
    alphaloc=linearptsFelas(i);
    peelloc=peel_f(i);
    if peelloc<alphaloc
        R_I_II=(anglestheo(i)+anglestheo(i-1))/2;
        break
    end
end
for i=1:length(anglestheo)
    peelloc=peel_f(i);
    eqUloc=intersectUFelas(i);
    if anglestheo(i)>R_I_II && peelloc>eqUloc
        R_II_III=(anglestheo(i)+anglestheo(i-1))/2;
        break
    end
end

%Determine Regions & Plot Data Points
R1loads=[];
R2loads_l=[];
R2loads_n=[];
R3loads_l=[];
R3loads_n=[];
R1angles=[];
R2angles_l=[];
R2angles_n=[];
R3angles_l=[];
R3angles_n=[];
nlocked=0;
nnestled=0;
% distinguish data by conformation and failure regions
%optional: scatter all non-prebuckled data 
for c=1:length(dat(:,1))
    phi=dat(c,1);
    if phi<21.9
        %scatter(dat(c,1),dat(c,2),30,'filled','MarkerFaceAlpha',1/3,'MarkerFaceColor',[0 0 0]...
        %    ,'MarkerEdgeColor',[0 0 0],'MarkerEdgeAlpha',2/3);
        R1loads=[R1loads dat(c,2)];
        R1angles=[R1angles dat(c,1)];
        nlocked=nlocked+1;
    elseif phi<=31.8 && dat(c,2)>=0.0256
        %scatter(dat(c,1),dat(c,2),30,'filled','MarkerFaceAlpha',1/3,'MarkerFaceColor',[255/255 112/255 101/255]...
        %    ,'MarkerEdgeColor',[0 0 0],'MarkerEdgeAlpha',2/3);
        R2loads_l=[R2loads_l dat(c,2)];
        R2angles_l=[R2angles_l dat(c,1)];
        nlocked=nlocked+1;
    elseif phi<=31.8 && dat(c,2)<0.0256
        %scatter(dat(c,1),dat(c,2),30,'filled','MarkerFaceAlpha',1/3,'MarkerFaceColor',[0/255 190/255 30/255]...
        %    ,'MarkerEdgeColor',[0 0 0],'MarkerEdgeAlpha',2/3);
        R2loads_n=[R2loads_n dat(c,2)];
        R2angles_n=[R2angles_n dat(c,1)];
        nnestled=nnestled+1;
    else
        %scatter(dat(c,1),dat(c,2),30,'filled','MarkerFaceAlpha',1/3,'MarkerFaceColor',[0/255 190/255 30/255]...
        %    ,'MarkerEdgeColor',[0/255 100/255 0/255],'MarkerEdgeAlpha',2/3);
        R3loads_n=[R3loads_n dat(c,2)];
        R3angles_n=[R3angles_n dat(c,1)];
        nnestled=nnestled+1;
    end
end
%scatter(elasticmaxdat(1,:)',elasticmaxdat(2,:)',30,'filled','MarkerFaceAlpha',128/255,'MarkerFaceColor',[255/255 112/255 101/255]...
%    ,'MarkerEdgeColor',[180/255 0/255 0/255],'MarkerEdgeAlpha',128/255);
R3loads_l=[R3loads_l elasticmaxdat(2,:)];
R3angles_l=[R3angles_l elasticmaxdat(1,:)];

%scatter(peelmaxdat(1,:)',peelmaxdat(2,:)',30,'filled','MarkerFaceAlpha',0/3,'MarkerFaceColor',[180/255 0/255 0/255]...
%    ,'MarkerEdgeColor',[70/255 0/255 130/255],'MarkerEdgeAlpha',3/3);
R3loads_n=[R3loads_n peelmaxdat(2,:)];
R3angles_n=[R3angles_n peelmaxdat(1,:)];

%determine averages for non-prebuckled data
nlocked=nlocked+length(elasticmaxdat);
percentlocked=nlocked/(nlocked+nnestled);
R3avgload=mean([R3loads_l R3loads_n]);
R2avgload=mean([R2loads_l R2loads_n]);
R23avgload=mean([R3loads_l R3loads_n R2loads_l R2loads_n]);

%Determine Mean and Range for given sections of the peeled data to draw the 
%Mean and Range lines and regions respectively
rangedat=sortrows([dat(:,1)' elasticmaxdat(1,:); dat(:,2)' elasticmaxdat(2,:)]',1);

binlims=[linspace(min(rangedat(:,1)),20,5) linspace(22,30,4) linspace(33,max(rangedat(:,1)),9)];
numbin=length(binlims);
datloc=zeros(6,numbin-1);
for i=1:(numbin-1)
    minlim=binlims(i);
    maxlim=binlims(i+1);
    if i==1
        datloc(1,i)=minlim;
    elseif i==(numbin-1)
        datloc(1,i)=maxlim;
    else
        datloc(1,i)=(maxlim+minlim)/2;
    end
    bin=[];
    for j=1:length(rangedat(:,1))
        angle=rangedat(j,1);
        if angle>=minlim && angle<=maxlim
            bin=[bin rangedat(j,2)];
        end
    end
    datloc(2,i)=mean(bin);
    datloc(3,i)=std(bin);
    if length(bin)==1
        datloc(4,i)=mean(bin);
        datloc(5,i)=mean(bin);
    else
        datloc(4,i)=min(bin);%mean(bin)-std(bin);%
        datloc(5,i)=max(bin);%mean(bin)+std(bin);%
    end
    datloc(6,i)=length(bin);
end

%Plot Average for All non-pre-buckled Data Points
plot(datloc(1,:),datloc(2,:),'Marker','none','Color',[130/255 130/255 130/255],'LineStyle','-','LineWidth',3);

%optional: plot range for all non-pre-buckled data
%jbfill(datloc(1,:),datloc(4,:),datloc(5,:),[200/255 200/255 200/255],[1 1 1],[0 0 0],119/255);
%plot(datloc(1,:),datloc(4,:),'Marker','none','Color',[130/255 130/255 130/255],'LineStyle','-');
%plot(datloc(1,:),datloc(5,:),'Marker','none','Color',[130/255 130/255 130/255],'LineStyle','-');


%Plot Peeling Model
plot(phiphi,PeelF(phiphi)./areacalculation(w,phiphi),'Marker','none','Color',[127/255 0 255/255],'LineStyle','-','LineWidth',1.5);

%Plot Tensile Strength of Tape
hline=ones(1,length(phiphi)).*TS_tape;
plot(phiphi,hline./areacalculation(w,phiphi),'Color','k','LineStyle','--');

%Determine & Plot Lap Shear Strength
LSSdata=importdata('LSS Raw Data.txt');
Data=LSSdata.data;
JL=Data(:,1)';
PeakForce=Data(:,2)';
JunctionArea=(JL).*w;
ff=@(m,x) m.*x;
[linreg]=fit(JunctionArea',PeakForce','poly1');
lapshearstrength=linreg.p1;
LSF=ones(1,length(phiphi)).*lapshearstrength;
phiphi=linspace(10,90,10000);
plot(phiphi,LSF,'k');

%Plot where Elastic Energy=Peeling Energy
plot(anglestheo,intersectUFelas,'Color',[255/255 193/255 37/255],'LineWidth',1.5);

%Plot where slope of elastic force=modulus of single tape (alpha)
plot(anglestheo,linearptsFelas,'Color',[255/255 126/255 24/255],'LineWidth',1.5);

%Plot Region Boundaries
line1_2=ones(1,1000).*R_I_II;
line2_3=ones(1,1000).*R_II_III;
yy=linspace(0,0.1,1000);
plot(line1_2,yy,'b','LineStyle','--');
plot(line2_3,yy,'b','LineStyle','--');

%Formatting the plot
xlabel('Junction Angle (\circ)','FontSize',14); 
xlim([10 90]);
ylabel('Max Load/Junction Area (MPa)','FontSize',14);
ylim([0 0.1]);
ax = gca;
ax.XAxis.MinorTick = 'on';
ax.XAxis.MinorTickValues = ax.XAxis.Limits(1):2.5:ax.XAxis.Limits(2);
ax.FontSize=14;
%%
figure;
prebuckle=categorical({'R_1', 'R_2','control'});
values=[percentlocked percentlocked4 percentlocked2].*100;
bar(prebuckle,values);